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ABSTRACT 

In  this  paper  we  discuss  a  numerical  approach  for  the  treatment  of  optimal  shape 
problems  governed  by  the  Euler  equations.  In  particular,  we  focus  on  flows  with  embedded 
shocks.  We  consider  a  very  simple  problem:  the  design  of  «i  quasi-one-dimensional  Laval 
nozzle.  We  introduce  a  cost  function  and  a  set  of  Lagrange  multipliers  to  achieve  the 
minimum.  The  nature  of  the  resulting  costate  equations  is  discussed.  A  theoretical  difficulty 
that  arises  for  cases  with  embedded  shocks  is  pointed  out  and  solved.  Finally,  some  results 
are  given  to  illustrate  the  effectiveness  of  the  method. 
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1  Introduction 


The  physical  structure  of  the  complex  flows  that  occur  in  aerodynamic  design  can  be 
predicted  by  reliable  numerical  simulations.  On  the  other  hand,  the  increasing  capability 
of  computers  to  perform  even  larger  calculations  radically  changes  the  aerodynamic  design 
process.  Indeed,  for  engineering  purposes,  if  one  can  predict  performance,  it  is  fundamental 
to  know  which  modification  of  an  aerodynamic  configuration  improves  performance. 

This  question  has,  of  course,  been  addressed  long  before  the  advent  of  computers  which 
has  led  to  a  broad  category  of  methods  known  as  inverse  design.  An  exhaustive  account 
of  the  historical  development  of  these  approaches  is  given  in  [1].  Here  it  is  suffice  to  say 
that  these  methods,  pioneered  by  Lighthill  [6],  require  knowledge  of  a  desirable  pressure 
or  velocity  distribution.  The  adequacy  of  the  distribution  chosen  is  dependent  on  the 
experience  of  the  designer;  the  resulting  shape  strongly  depends  on  this  choice.  An  original 
example  of  such  an  approach  is  found  in  [3] 

The  numerical  approach  that  we  will  use  in  this  paper,  lift  the  dependence  on  heuristic 
choices  of  the  desirable  distribution,  allowing  the  imposition  of  constrains  to  be  satisfied 
by  the  solution  found.  The  numerical  simulation  of  the  flow  and  a  numerical  optimization 
code  are  coupled.  The  optimization  code  calculates  the  best  perturbation  to  the  geometry 
to  decrease  a  cost  function.  The  geometry  itself  is  described  by  a  set  of  shape  coefficients. 

The  optimization  code  can  be  devised  in  one  of  several  ways.  A  common  approach  is 
to  perturb  one  shape  coefficient  at  a  time  and  compute  the  derivative  of  the  cost  function 
with  respect  to  this  coefficient  as  a  finite  difference.  Although  such  codes  are  simple  to 
devise,  the  procedure  is  costly  and  can  introduce  large  errors.  In  a  further  evolution  of  this 
approach,  an  equation  is  first  calculated  for  the  derivative  of  the  cost  function  with  respect 
to  the  shape  coefficient  and  then  solved  numerically.  An  equation  must  be  solved  for  each 
shape  coefficient.  A  recent  application  of  this  method  to  a  two-dimensional  supersonic 
problem  is  found  in  [2]. 

The  approach  presented  in  this  paper  is  a  classical  optimal  control  method.  We  will 
introduce  costate  variables  (Lagrange  multiplier)  to  achieve  a  minimum.  This  method  has 
been  successfully  applied  in  the  design  of  an  airfoil  in  a  subsonic  potential  flow  [10]. 

Here,  we  consider  a  flow  with  embedded  shocks  where  the  governing  equations  are  the 
Euler  equations.  We  show  how  to  derive  an  analytical  expression  of  the  cost  function 
derivatives  with  respect  to  the  shape  coefficients.  For  this  purpose,  we  solve  only  one  set 
of  costate  equations.  In  [9],  some  difficulties  are  outlined,  and  a  method  to  avoid  them 
is  proposed.  In  the  present  approach  an  optimal  shape  can  be  found  for  problems  with 


embedded  shocks,  without  additional  complications.  A  careful  examination  of  the  structure 
of  the  costate  equations  suggests  a  method  for  integrating  them  with  a  robust  algorithm 
developed  for  fluid  dynamics  purposes. 

A  comparison  of  optimization- based  approaches  for  aerodynamics  design  problems  is 
given  in  [7],  although  some  results  for  flows  with  embedded  shocks  have  been  questioned 
recently  in  [9], 
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2  Problem  Statement 


We  investigate  a  new  method  for  aerodynamic  design  and  optimization  based  on  the  Eu¬ 
ler  equations.  In  order  to  demonstrate  the  ideas  of  the  method,  a  very  simple  problem 
is  considered:  the  design  of  a  Laval  nozzle  assuming  inviscid,  quasi-one-dimensional  flow. 
The  optimization  problem  consists  of  finding  a  set  of  design  variables  (in  this  case  shape 
parameters)  that  minimize  some  cost  function,  e.g.,  a  desired  pressure  or  velocity  distribu¬ 
tion  along  the  centerline  and  possibly  requiring  that  some  side  constraints  be  satisfied.  We 
attack  the  optimization  problem  with  the  adjoint  method.  The  adjoint  method  introduces 
a  new  set  of  equations  and  unknowns  that  are  solved  together  with  the  flow-field  equations. 
To  better  understand  how  to  obtain  a  solution  of  the  adjoint  equations,  the  properties  of 
these  equations  are  discussed. 

Let  the  extension  of  the  Laval  nozzle  in  physical  space  be  Q  =  [0,/].  An  energylike 
functional  denoted  by  £,  is  the  cost  function  we  want  to  satisfy.  An  optimal  shape  of  the 
nozzle  is  reached  when  we  meet  the  necessary  conditions  for  a  minimum  of  £.  Let 

E  =  \Jo^P~P^2dX 

where  p*  is  a  target  pressure  distribution  along  the  abscissa  x  and  p  is  the  pressure  field 
for  the  present  geometry.  The  choice  of  the  functional  does  not  affect  the  generality  of  the 
method  once  the  dependence  of  £  with  respect  to  the  flow-field  variables  is  determined. 
Nevertheless,  the  choice  of  the  functional  itself  (e.g.,  |p  —  p*|  instead  of  (p  — p*)2)  can  affect 
the  performance  of  the  optimization  algorithm  by  changing  the  curvature  of  the  energylike 
surface. 

In  the  case  of  a  quasi-one-dimensional  flow,  if  we  denote 

p  =  density 
u  =  velocity 
e  =  specific  total  energy 
p  =  pressure 
a  =  speed  of  sound 
h  =  height  of  the  channel 
7  =  specific  heat  ratio 


then  the  Euler  equations  for  unsteady  flows  reduce  to 

U,  +  F,  +  Q  =  0  (1) 

where 

(  P  \  (  P*  \  (  P«P  \ 

U=  pu  F  =  I  p  +  pu2  Q  =  p«2/? 

\  pe  /  \  u(pe  +  p)  /  \  u(pe  +  p)0  / 
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/3  =  A*/A,  and  p  =  np(2e  —  u2).  In  the  following  derivation,  we  use  the  homogeneity 


property 

F  =  A(U)U 

(2) 

where 

dF  ,  s 

• 

au=A(u) 

(3) 

and 

/ 
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A(U)  = 
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The  source  term  Q  can  be  written  to  display  its  dependence  on  U;  in  fact,  the  multi¬ 
plication  can  be  carried  out  to  show  that 


Q  =  S(U)U 

where 


(4) 


S(U)  =  — 
{  }  dV 


0  1  0  \ 
-u2  2u  0 
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The  first  and  the  third  tow  of  A(U)  are  proportional  to  the  first  and  third  row  of  S(U), 
respectively  .  Furthermore  because  d(pu2)  =  —  u2  -dp+2u-d(pu)  and  pu2  =  —u2  p  +  2u  pu, 
it  follows  that  S(U)  =  dQ/d\J. 

We  refer  to  the  solution  of  the  above  equations  as  the  analysis  problem.  The  bound¬ 
ary  conditions  for  this  problem  must  be  chosen  such  that  the  problem  is  well  posed.  In 
particular,  we  will  consider  the  inlet  flow  with  a  constant  total  pressure  and  entropy,  i.e., 
Pin  =  P«n(  1  +  K-M2nyt'y~x  =  Constant,  Pin/ pin  =  Constant.  At  the  outlet,  if  the  flow  is 
subsonic,  the  static  pressure  is  fixed  as  p ^  =  Constant. 


3  Adjoint  Formulation  for  a  Shockless  Case 

The  design  problem  can  be  thought  of  as  a  search  for  a  minimum  of  a  functional  under 
constrains.  Let 

£(<*,)  =  £  +  /  AT(Fr  +  Q )dx  (6) 

where  Ar  is  an  arbitrary  vector  with  components  (Aj(x),  Aj(a;),  A3(z)),  fi  is  the  domain 
[0, /],  and  are  shape  coefficients  that  define  the  geometry  of  the  nozzle,  for  example  by 
h(a,,x)  =  J2iaifi{x)  with  /,(x)  a  generic  function  of  x.  Since  in  the  steady  state  the 
Euler  equations  must  be  satisfied  everywhere  in  the  domain,  the  functionals  £  and  £  are 
identical.  If  we  increase  the  shape  coefficients  by  eai,  then  the  latter  functional  will  change 
by  an  amount,  say,  sSC.  The  other  quantities  will_change  in  the  same  way:  0  ^  fl  +  efl, 
p  =»  p  +  ep,  and  U(x)  =>•  U(x)  ■+•  eU(x),  where  U  =  (p,pu,pe)T.  By  using  eqs.(3)  and 
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(5),  we  obtain  also  F(x)  =>  F(x)  +  eA(U)U(x)  and  Q(x)  =»  Q(x)  +  eS(U)U(x).  If  we 
substitute  these  relations  into  eq.(6)  and  retain  only  the  first-order  terms,  we  obtain 

SC  =  [‘ftp  -  p')dx  +  /'  At{[A(U)U]x  +  S(U)U}dx  +  /'  ^HlHdx  (7) 
Jo  Jo  Jo  P 

In  the  above  equation,  =  (hxh  —  hhx)/h2.  With  the  notation  Cj(x)  =  ( hfix  —  hxfi)/h 2, 
the  last  integral  of  eq.(7)  can  be  written  as 

in  which  the  substitution  ft  =  5^  is  made. 

Let  us  integrate  by  parts  the  second  integral  in  eq.(7),  with  A  =  A(U)  and  S  =  S(U) 
We  obtain 

At[(AU)i  +  SUjdx  =  [ATAU]{,  -  [‘  AjAUdx  +  [‘  ArS Udx  (9) 

Jo  Jo  Jo 

The  first  term  on  the  right-hand  side  of  the  above  equation  drops  for  a  suitable  choice  of 
A  at  the  boundaries.  The  boundary  conditions  for  U  are  complementary  to  those  for  A, 
in  the  sense  that  ATAU  =  0  yields  a  homogeneous  linear  system  in  A  whose  rank  depends 
on  the  number  of  boundary  conditions  for  U. 

For  this  test  case,  at  the  inlet  we  have 

p°n  =  Constant 

which  implies  that 

p  =  —puu 

If  the  entropy  is  fixed,  then  the  specific  total  enthalpy  is  also  constant;  therefore, 

7 — — I-  ^u2  =  Constant 
(7  -  1  )p  2 

We  conclude  that 

[AU]0  =  [pu,  upu,  +  2U^^T 

Hence,  the  suitable  choice  at  the  inlet  for  A  is 

Aj  +  UA2  +  7  +  - u 2  A3  =  0  (10) 

1(7-1)^  2 

At  the  outlet  for  a  subsonic  case  with  a  given  pout ,  we  obtain 


(AU]i  =  {pu,2tipu  —  u2p, 


- 's - (-  —u‘ 

(7  ~  1  )P  2 


3  ,1  _ 


pu  —  u 


(7  - 1  )p 


4 


which  leads  to 


Ai  ■+■  2uAa  -f 
t*A3  + 


IP 


_  ,  3  2 

(7  ~  1  )P  2U  J 
IP 


+  u‘ 


A3  =  0 

a3  =  0 


(11) 


1(7 -l)p 

If  the  outflow  is  supersonic,  then  no  boundary  conditions  are  required  for  U;  therefore, 

A  =  0  (12) 


identically  at  the  outlet. 

If  we  take  boundary  conditions  (10)  and  (11)  and  eq.(9)  into  account,  eq.(7)  can  be 
rewritten  as 


SC 


=  J UT|-  AtA.  +  StA  +  | £(p  -  p-)]dz  + 


E  */ 

i  Jo 


c,AtSU 

p 


dx 


(13) 


where  dp/d U  =  ^^(u2,  —  2u,  2).  If  we  select  A  such  that 

-ATA.  +  SrA  +  ^(P-p-)  =  0 


(14) 


eq.(13)  becomes 


in  which  we  recognize 


,  CjArSU 


P 


dx 


dC  _  r‘  a  ATSU 

doti  Jo  P 


(15) 


Suppose  that  the  flow  field  is  known,  such  that  all  the  variables,  dependent  upon  U,  are 
fixed.  If  we  solve  eq.(14)  with  the  appropriate  boundary  conditions  for  A  and  substitute 
the  results  in  eq.(15),  then  we  obtain  a  formulation  for  the  gradient  of  the  Lagrangian.  The 
problem  then  is  reduced  to  finding  the  solution  of  a  linear  system  of  equations  in  A  with 
homogeneous  boundary  conditions.  Note  that  A  =  0  is  a  solution  if  p  =  p*,  which  is  a 
sufficient  condition  for  SC  =  0.  If  at  the  minimum  p  ^  p*  although  the  integral  in  eq.(15)  is 
0,  then  in  general  A  ^  0.  The  discussion  thus  far  leads  to  some  basic  questions  about  the 
well  posedness  of  the  eq.(14)  with  the  boundary  conditions  given  by  eqs.(10),  (11)  or  (12) 
and  about  the  existence  and  uniqueness  of  the  solution.  To  actually  determine  a  solution 
of  this  system,  examine  eq.(14),  embedded  in  time  as 


±  At 


ATA*  +  StA  + 


dp_ 

du 


(p-p*)  =  0 


(16) 


in  which  we  must  choose  for  the  proper  sign  for  the  time  derivative. 
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The  matrix  AT  has  the  familiar  eigenvalues  u  —  a,u  and  u-fa.  At  the  boundaries,  if  we 
choose  the  positive  sign  in  eq.(16),  for  each  boundary  condition  there  will  be  an  incoming 
characteristic,  such  that  the  problem  is  well  posed.  Another  motivation  for  this  choice 
is  that  if  we  add  to  eq.(l)  the  diffusive  term,  when  integrated  by  parts  twice,  it  would 
append  to  eq.(14)  a  second  derivative  term  in  x  with  the  same  sign  this  term  had  in  the 
flow  equations.  Therefore,  if  the  negative  sign  in  eq.(16)  is  adopted,  then  an  ill-posed  heat 
equation  would  result.  In  conclusion,  note  that  the  costate  equation 

A<  -  AtAz  +  SrA  +  ^(p  -  p*)  =  0  (17) 

has  an  upside-down  characteristic  pattern  with  respect  to  the  time-dependent  flow  equation. 
If  we  consider  a  transonic  nozzle,  in  the  throat  area,  the  eigenvalue  u  —  a  goes  continuously 
through  zero;  the  flow  undergoes  an  expansion  through  a  transonic  fan.  On  the  other  hand, 
the  behavior  of  the  same  family  of  characteristics  for  eq.(17)  shows  a  shocklike  pattern.  The 
two  characteristic  patterns  are  illustrated  in  figs.  1(a)  and  1(b). 

For  a  reason  that  will  be  clear  later,  some  ambiguous  jump  conditions  for  this  shock 
will  be  derived.  First  consider  a  simple  equation  of  the  kind  $t  —  x$x  =  0,  with  0  =  $(x), 
x  €  [—1,1],  and  boundary  conditions  $(— 1)  =  $/,  $(1)  =  $r.  The  characteristics  at  the 
boundaries  show  that  this  problem  is  well  posed  and  independent  of  the  initial  conditions, 
for  a  large  t ,  the  solution  will  be  a  step  function  $(x)  =  for  x  €  [— 1,0[  and  4>(x)  =  $r 
for  x  €]0, 1].  Note  that  the  jump  at  zero  is  solely  determined  by  the  boundary  conditions 
and  that  the  steady  solution  will  be  reached  for  t  — ►  oo  because  of  the  nearly  vertical 
characteristics  next  to  x  =  0. 

The  structure  of  the  solution  to  this  equation  is  similar  to  that  underlying  eq.(17)  for 
which  the  analysis  hides  somehow  this  behavior.  A  solution  to  eq.(14)  can  be  written  in 
the  form 

A(x)  =  C(x)  +  [A ]H(x  -  xth) 

where  C(x)  €  C1,  [A]  is  the  jump  at  the  shock,  H(x  —  xth)  is  the  Heaviside  function,  and  xth 
is  the  location  of  the  shock  for  A  (i.e.,  the  throat  of  the  nozzle).  If  we  substitute  in  eq.(14), 
because  6(x)  (which  is  the  Dirac  measure  of  x)  is  the  derivative  of  H(x)  with  respect  to  x, 
we  obtain  at  x  =  xtj, 

AT[A]6(x  xt*,)  =  0 

in  which  all  the  negligible  terms  have  been  dropped.  Note  that  the  presence  of  source  terms 
in  eq.(14)  does  not  affect  the  derivation  of  the  jump  conditions. 

Finally,  because  at  x  =  xth  det  A  =  0,  we  obtain  a  nontrivial  solution  for  the  system 

Ar(A]  =  0  (18) 

which  yields  only  two  jump  conditions;  the  third  jump  condition  depends  on  the  boundary 
conditions.  In  this  sense  this  problem  has  ambiguous  jump  conditions. 

If  the  only  solution  of  the  homogeneous  problem  associated  with  eq.(14),  i.e. 

-  ATAX  +  STA  =  0  (19) 
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and  its  boundary  conditions,  is  A  =  0,  then  thfc  solution  with  a  nonhomogeneous  source 
term  is  unique,  since  eq.(14)  is  linear.  Let  us  consider  a  subsonic  case,  with  det  A  ^  0 
everywhere  in  the  domain.  A  general  solution  of  eq.(19)  can  be  written  as 

A  =  Ao  *£<*'>-*'* 

and,  together  with  boundary  conditions  (10)  and  (11),  this  implies  A  =  0  on  the  domain. 
In  the  transonic  case,  since  detA  =  0  at  x  =  x/^,  we  split  the  problem  into  two  domains, 
such  that  in  the  subsets  [0,Xt/,[  and  ]xt^,  /],  detA  ^  0.  The  solution  will  be 

A,„6  =  A'0efo,h(AT)-,sTjx  for  x  €  [0,xtA( 

and 

*  *//  (AT)-‘SrAr  ,  ,  „ 

Asuper  =  A0eJ**»  for  x  eJxtA,  /]. 

Again,  if  we  account  for  the  inlet  boundary  conditions  (10),  the  jump  conditions  from 
eq.(18),  and  the  outlet  boundary  conditions  (12),  then  the  solution  is  A  =  0. 

In  summary,  we  have  derived  an  analytic  formulation  for  the  gradient  of  the  Lagrangian 
in  eq.(6)  with  respect  to  the  geometry.  Furthermore,  we  have  shown  that  this  representation 
is  unique  in  the  sense  discussed  above. 


4  Costate  Equations  for  a  Shock  Case 


Until  now,  we  have  limited  our  investigation  to  shockless  nozzles  to  avoid  certain  difficulties 
that  we  will  discuss  here.  One  problem  is  that  eq.(l)  and,  therefore,  eq.(16)  are  not  defined 
at  the  shock.  This  problem  is  overcome  by  extending  the  solution  space  of  U(x)  to  a  set  of 
generalized  functions,  such  that  eq.(l)  will  reduce  to  the  Rankine-Hugoniot  jumps  at  the 
shock.  A  more  subtle  shortcoming  is  better  understood  with  the  aid  the  following  example. 
Consider  a  simple  equation  of  the  kind  $t  +  (//(x)  —  l/2)$r  =  0  that  is  defined,  for  example, 
on  ft  =  [—1,1].  The  characteristics  pattern  (fig.2),  shows  the  necessity  of  some  boundary 
conditions  on  both  sides  of  the  the  discontinuity  to  ensure  the  existence  of  a  steady-state 
solution,  regardless  of  the  boundary  conditions  at  the  ends  of  the  domain  ft. 

Now,  eq.(17)  can  be  rewritten  to  reflect  its  characteristics  pattern 


where 


dp 


P i  At  -  DPTAr  +  P  S  A  +  P T-^~{p  -  p*)  =  0 


P  = 


1 

u  —  a 


u 

\  (e  +  p)  -  ua  ±u2  (e  +  p)  +  ua 


1 

u  +  a 


(20) 


is  the  matrix  of  the  right  eigenvectors  of  A(U),  and 
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/  u-a  0  0  \ 

D  =  0  u  0 

\  0  0  u+a ) 


At  the  shock  wave,  the  characteristic  that  corresponds  to  the  eigenvalue  u  —  a  undergoes 
a  jump  in  speed  of  the  kind  described  in  the  example.  If  this  point  is  considered  inside 
the  domain  of  calculation,  then  some  condition  is  needed  to  update  the  solution  in  time. 
A  boundary  condition  is  needed  at  this  point  to  continue  the  calculation  to  the  left  of  the 
shock  wave.  We  cannot  try  to  derive  some  boundary  conditions  for  this  point,  as  explained 
for  the  inlet  and  the  outlet.  No  speculation  about  the  perturbation  U  at  the  shock  is 
possible  because  a  perturbation  that  is  solely  dependent  on  the  shape  coefficients  a,  and 
the  flow  equations  would  be  chosen  arbitrarily.  No  other  constrains  exist.  For  example, 
the  application  of  the  Rankine-Hugoniot  jumps  to  U  on  each  side  of  the  shock  would  be 
equivalent  to  assuming  that  the  shock  does  not  change  position  regardless  of  the  value  of 
The  integrals  in  eq.(7)  are  split  in  two,  and  the  integration  is  carried  between  [0, 
and  ]x3h,  /]  as 

SjC  —  6C,\  -f- 


where 


SC,  =  [ArAUl0'-'  +  £"  Ur|-ATA,  +  StA  +  -  ,>•)]* 

__  r.K  c,AtSU 

?QiX 


+ 


0 


-dx 


(21) 


and 


dp 


sc ,  =  [ATAUt„+jf_U7'(-ATAI  +  STA+^(p-p-))<ii 


+ 


X,h 

_  r'  c,ATSU 


!>/ 


dx 


(22) 


A  suitable  choice  for  the  Lagrange  multiplier  is  to  take  A  =  0  on  both  sides  ofthe  shock 
such  that  A  is  continuous.  This  selection  frees  us  from  imposing  a  condition  on  U,  because 
the  addendum  [ATAU]  in  eqs.(21)  and  (22)  drops  anyway  at  the  shock.  At  the  shock,  we 
will  have  three  characteristics  that  deliver  the  information  A  =  0  to  the  left  domain,  and 
one  that  delivers  it  to  the  right. 

We  have  examined  also  another  possible  interpretation  of  eqs.(21)  and  (22).  Assume 
that  for  some  perturbation  of  the  shape  coefficients  or;,  the  shock  properties  do  not  change, 
i.e.,  the  jump  in  total  pressure  across  it,  is  essentially  constant,  which  will  be  the  case 
for  a  weak  shock.  If  we  assume  that  the  total  enthalpy  is  constant,  then  the  field  to 
pa  formulationthe  right  of  the  shock  can  be  regarded  as  a  subsonic  nozzle  governed  by 
eq.(17)  with  boundary  conditions  (10)  and  (11).  To  the  left  of  the  shock,  the  flow  behaves 
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as  a  supersonic  outlet,  where  we  must  impose  the  proper  boundary  condition,  eq.(12),  as 
discussed  earlier.  With  this  approach,  8C\  =  0  and  8C2  =  0  independently  at  the  minimum. 

The  two  approaches  presented  to  handle  the  shock  are  not  significantly  different.  In 
the  numerical  calculations  that  follow,  we  have  implemented  the  second  set  of  boundary 
conditions. 


5  Numerical  Approach 


The  flow-field  solution  is  obtained  by  introducing  a  discrete  grid  ( xn,t* )  =  (xo  +  nAx,t0  + 

At*),  where  Ax  is  constant  and  At*  changes  to  satisfy  the  CFL  condition.  The  con¬ 
servative  variables  U(x)  are  computed  at  the  cell  centers  and  integrated  in  time  with  a 
three-stage  Runge-Kutta  scheme,  as  explained  in  [5].  In  this  implementation,  we  interpo¬ 
late  U(x)  to  the  cell  faces  by  using  characteristic  differences  and  a  minmod  limiter.  The 
flux  derivative  in  eq.(l)  is  then  computed  using  an  approximate  Riemann  solver.  See  [4], 
Away  from  discontinuities,  the  scheme  is  second  order  accurate. 

Depending  on  the  case  considered,  the  solutions  of  the  costate  eq.(14),  are  sought  as 
steady  results  of  eq.(17),  with  boundary  conditions  applied  as  explained  in  the  previous 
section.  Although  eq.(17)  is  linear,  it  presents  some  numerical  difficulties  because  of  the 
characteristics  pattern  at  the  throat  and  at  the  shock  location.  In  particular,  consider 
eq.(20)  discretized  over  the  same  uniform  grid  of  the  analysis  problem  with  spacing  Ax.  If 
we  denote  by  A(-)u  the  finite  increments  of  the  function  (•)  with  respect  to  the  superscripted 
variable,  we  have,  at  each  grid  point 

pT^-DpT^7  +  pTsTA  +  pT^(p-',')  =  0  (23) 

Define  the  local  increment  AW  =  PTAA;  with  this  notation  eq.(23)  becomes 


AW' 


-D 


AW1 


+  PTSTA  +  PTJ^(p-p-)=0. 


(24) 


At  Ax 

This  equation  describes  the  signals  that  propagate  along  the  characteristics;  therefore, 
the  increment  AW1,  is  one-sided  depending  on  the  sign  of  the  corresponding  propagation 
speed.  Note  that  for  this  equation  it  would  be  impossible  to  use  a  conservative  scheme 
since  no  conservation  law  exists  to  satisfy.  The  integration  in  time  is  made  by  explicit  time 
stepping.  The  scheme  is  first-order  accurate. 

At  the  throat  of  the  nozzle,  u  —  a  =  0  and  =  0;  hence,  the  first  row  of  eq.(24)  reduces 
to 


AW1 

At 


(p-p*)  =  0 


The  eigenvalue  u  —  a  has  been  shown  to  vanish  at  the  throat.  For  a  grid  point  in  the 
neighborhood  of  the  throat,  this  singularity  can  lead  to  unbounded  grows  for  Lambda, 
depending  on  the  nature  of  the  source  term  in  the  above  equation.  Because  the  other  two 
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characteristics  are  nonzero  at  the  throat,  this  error  can  degrade  the  entire  calculation.  To 
avoid  this  problem,  we  approximate  u  —  a  with  its  value  at  the  neighboring  point  on  the 
side  from  which  the  vanishing  characteristic  propagates. 

Since  with  a  Godunov-type  solver  the  shock  is  resolved  with  three  grid  points,  we  must 
decide  on  which  of  these  to  impose  the  boundary  conditions  for  A.  We  must  consider  that 
the  middle  point  of  the  three  cells  on  which  the  shock  is  solved,  is  almost  sonic;  if  the 
boundary  conditions  for  A  at  a  subsonic  inlet  (eq.(10))  were  imposed  at  this  grid  point, 
then  the  convergence  rate  to  the  steady  solution  would  be  considerably  slower.  For  this 
reason  we  impose  the  condition  for  supersonic  outlet  A  =  0,  on  the  middle  grid  point, 
eliminating  it  from  of  the  computation. 

Another  remark  should  be  made  in  regard  to  the  order  of  magnitude  of  the  residu¬ 
als  of  eq.(23),  for  which  we  can  consider  the  solution  steady.  Close  to  the  minimum,  the 
gradients  in  eq.(l5)  are  almost  zero;  nevertheless  the  optimization  algorithm  requires  a 
careful  computation  of  these  values,  such  that  in  order  to  consider  the  time-dependent  so¬ 
lution  converged,  the  residuals  must  be  some  orders  of  magnitude  smaller  than  the  gradient 
components. 

In  the  results  that  follow,  we  used  a  representation  of  the  nozzle  geometry  defined  by 
h  =  a^X  +  a-i/X  +  0*3,  where  X  =  x  -f  10-3;  this  representation  allows  two  independent 
design  variables  because  (3  =  hx/h. 

In  this  work,  we  do  not  address  the  methods  to  accelerate  the  numerical  scheme  to 
obtain  an  optimal  shape;  the  strategy  used  to  achieve  the  minimum  of  the  functional  is 
straightforward: 

1.  Start  with  a  first  guess  for  the  shape  coefficients. 

2.  Solve  the  flow  equation. 

3.  Solve  the  costate  equation  with  the  values  computed  in  step  2  for  the  flow  field. 

4.  Update  the  shape  coefficient  with  a  gradient-based  criterion. 

5.  Restart  the  procedure  from  step  2  until  the  gradient  is  zero. 

To  update  the  shape  coefficients  a  BFGS  algorithm  was  used.  See  [8].  In  some  cases, 
as  we  will  discuss,  we  used  an  inefficient,  but  robust,  algorithm  that  simply  makes  a  line 
search  for  a  zero  of  the  gradient. 


6  Discussion  of  the  Results  and  Conclusions 

The  values  for  the  functional  £,  computed  with  an  analytical  solution  of  eq.(l),  are  shown 
in  fig.(3).  The  numerical  values  of  the  analytical  solution  are  computed  on  the  same  grid 
presented  above;  then,  the  functional  is  computed  by  a  trapezoidal  approximation.  The 
discrete  functional,  which  is  a  result  of  the  trapezoidal  integration,  shows  some  disconti¬ 
nuities  and  a  local  minimum  that  disappears  as  the  number  of  grid  points  increases.  As 
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the  mesh  is  refined,  the  number  of  the  discontinuities  increases,  while  the  jumps  become 
smaller. 

Note  that  if  the  dependence  of  the  geometry  on  the  shape  coefficients  is  smooth,  then 
the  functional 

d£  [dp 

torLfoi{r-p)dx 


is  always  defined  with  the  assumption  that  dp/dati  is  defined  everywhere  except  at  a  finite 
number  of  points.  This  assumption  is  reasonable  because  the  solution  of  eq.(l)  changes 
smoothly  with  the  geometry.  For  this  reason,  even  if  no  certainty  exists  that  the  solution 
depends  monotonically  on  the  shape  coefficients,  this  behavior  can  be  the  interpreted  in 
this  way.  Suppose  that  the  integrand  of  the  functional  can  be  represented  by  a  simple 
rectangular  function.  If  in  attaining  the  minimum  the  area  under  the  curve  decreases  and  its 
“height”  increases,  the  functional  will  eventually  increase  before  the  edges  of  the  rectangle 
pass  another  grid  point,  because  the  mesh  resolution  is  not  sufficient.  The  functional  will 
exhibit  a  local  minimum  and  a  subsequent  discontinuity. 

A  method  that  derives  a  formulation  for  the  gradient  of  £  from  a  discrete  approximation 
of  the  functional  will  ohtain  meaningless  solutions  as  a  result  of  the  discontinuities  of  the 
discrete  functional,  such  that  no  optimization  algorithm  alone  could  anyway  get  to  the 
minimum.  In  the  present  formulation,  an  approximate  representation  of  the  analytical 
gradient  of  the  functional  is  derived.  For  this  reason,  the  approximation  of  the  analytical 
gradient  will  be,  at  most,  affected  by  discontinuities  due  to  the  discretization  and  will  be 
always  monotonic  (if  the  analytical  functional  does  not  change  curvature)  and  bounded.  See 
fig.4.  In  figs. 5  and  6,  we  present  two  sets  of  results  in  which  the  target  pressure  distribution 
was  generated  with  the  same  h(x)  that  was  used  in  the  optimization  procedure.  In  fig.5, 
the  target  pressure  distribution  is  obtained  starting  from  a  subsonic  first  guess.  This  result 
shows  the  effectiveness  of  the  method.  Fig.6  shows  that  we  can  achieve  the  optimum  from 
both  sides  of  the  discontinuity. 

In  general,  when  p*(x)  is  fixed,  the  minimum  of  £  is  reached  for  different  values  of  the 
shape  coefficient  a,.  These  different  values  depend  on  whether  one  considers  the  analytical 
or  the  discretized  functional,  even  if  the  results  are  converged  on  the  grid.  If  the  target 
solution  can  be  attained  exactly,  then  both  values  coincide.  The  gradient  calculated  after 
the  proposed  derivation,  will  still  depend  on  the  discretization  through  the  nonhomogeneous 
term  in  eq.(14).  In  fig.7,  we  show,  for  a  case  in  which  p*  cannot  be  reached  exactly,  that 
the  distance  between  the  two  minima  becomes  approximately  half  when  the  grid  resolution 
is  doubled.  This  result  supports  the  hypothesis  that  the  minimum  calculated  through  the 
analytical  gradient  will  indefinitely  approach  to  the  actual  minimum  as  the  grid  is  refined. 

In  conclusion,  a  method  has  been  presented  to  calculate  the  gradient  components  of  a 
generic  functional,  in  which  (regardless  of  the  number  of  the  shape  coefficients)  only  one 
linear  costate  equation  must  be  solved.  The  minimum  computed  in  this  way  differs  from 
the  minimum  of  the  discrete  functional;  however  these  minima  indefinitely  approach  as  the 
grid  is  refined. 
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Figure  1:  Characteristic  patterns,  (a)  Transonic  expansion,  (b)  Correspondent  shocklike 
structure  for  costate  equations. 


Figure  2:  Characteristic  pattern  for  +  [H(x)  -  1/2]#*  =  0.  Boundary  conditions  are 
needed  on  both  sides  of  discontinuity. 
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(a)  (b) 


Figure  3:  Analytical  solution  calculated  with  shape  function  h(x)  =  a(x3  —  x2)  +  105; 
functional  has  been  calculated  with  trapezoidal  approximation,  (a)  Solution  distributed 
over  20  grid  points,  (b)  Solution  distributed  over  80  grid  points. 
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(a)  (b) 


Figure  4:  Result  shown  is  obtained  with  h(x)  =  a\X  +  0.3/X  +  10.  For  each  value  of  aj 
f  :w  equation  is  solved  by  a  Godunov  like  scheme,  and  gradient  is  calculated  as  proposed 
in  this  paper.  Each  time  the  shock  goes  through  grid  point,  discretized  functional  does 
not  have  a  monotonic  derivative;  gradient  has  discontinuities,  but  is  still  monotonic,  (a) 
Functional,  (b)  Gradient. 
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Pressure 


Figure  5:  Target  solution  is  h(x)  =  2.5A'+0.3/A'+5.  First  guess  is  h(x)  =  2X +0.52/A"+5. 
Shape  function  for  whirl  minimum  is  soug!  \  is  h(x)  =  aj  X  +  Q2/X  +  5,  and  optimization 
algorithm  used  is  BFGS.  Starting  value  of  the  functional  is  of  order  10-3.  At  minimum  it 
is  of  order  10-9.  Gradient  components  are  of  order  10-12  at  minimum. 
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Figure  6:  Target  solution  is  h(x)  =  2.5A’+0.3/A’+10.  Firs,  guess  is  h(x)  =  5X+0.3/X+10. 
Shape  function  for  which  minimum  is  sought  is  h(x)  =  a\X  +  aijX  + 10,  and  optimization 
algorithm  used  is  BFGS.  Starting  value  of  functional  is  of  order  of  10-3.  At  minimum  it  is 
of  order  10~9.  Gradient  components  are  of  order  10-12  at  minimum. 
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